Variational cluster approach to ferromagnetism in infinite dimensions and in 

one-dimensional chains 
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The variational cluster approach (VGA) is applied to study spontaneous ferromagnetism in the 
Hubbard model at zero temperature. We discuss several technical improvements of the numerical 
implementation of the VGA which become necessary for studies of a ferromagnetically ordered 
phase, e.g. more accurate techniques to evaluate the variational ground-state energy, improved local 
as well as global algorithms to find stationary points, and different methods to locate the magnetic 
phase transition. Using the single-site VGA, i.e. the dynamical impurity approximation (DIA), the 
ferromagnetic phase diagram of the model in infinite dimensions is worked out. The results are 
compared with previous dynamical mean-field studies for benchmarking purposes. The DIA results 
provide a unified picture of ferromagnetism in the infinite-dimensional model by interlinking difi'erent 
parameter regimes that are governed by different mechanisms for ferromagnetic order. Using the 
DIA and the VGA, we then study ferromagnetism in one-dimensional Hubbard chains with nearest 
and next-nearest-neighbor hopping tg- In comparison with previous results from the density-matrix 
renormalization group, the phase diagram is mapped out as a function of the Hubbard-[/, the 
electron filling and i2. The stability of the ferromagnetic ground state against local and short-range 
non-local quantum fiuctuations is discussed. 

PACS numbers: 71.10.Fd, 71.10.Hf, 75.10.-b, 75.10.Jm 



I. INTRODUCTION 

Itinerant fcrromagnctisni of nanomcter-sizcd 
transition-metal systems deposited on non-magnetic 
surfaces has attracted much attention recently. It is 
a fascinating physical but also technological vision to 
control the geometrical arrangement of the nanosystem 
on the atomic scale while studying its magnetic prop- 
erties with atomic resolutionii This provokes new and 
exciting questions. It is highly interesting, for example, 
to understand how many atoms are necessary and how 
these atoms should be arranged geometrically to ensure 
a stable ferromagnetic state. 

One important issue for magnetic nanosystems is their 
stability against thermal fluctuations i^ This is mainly de- 
termined by anisotropics. Anisotropic contributions to 
the total energy of a magnetic system can be several or- 
ders of magnitude higher at a surface or for a small clus- 
ter or chain as compared to a three-dimensional bulk of 
the same material. Nevertheless, the anisotropy strength 
is usually still much smaller than the exchange coupling 
and can thus be safely disregarded for the question of 
whether or not a ferromagnetic ground state is existing. 

The stability of a ferromagnetic ground state is, there- 
fore, a matter of quantum fiuctuations. Opposed to 
antiferromagnetic order, for example, the order param- 
eter is a conserved quantity in the case of ferromag- 
netism. It is thus mainly the quantum fluctuations of 
the paramagnetic state which are important and which 
the spin-polarized state is competing with. 

The ground state of itinerant systems;^ such as mono- 
atomic chains of 3d transition nietals)^ is of particular 
interest from a theoretical point of view. Besides the 
subtle interplay between the kinetic energy of the itin- 
erant electrons and their Coulomb interaction, geomet- 



rical constraints come into play additionally. Even for a 
bulk system, however, and even for the most elementary 
models of itinerant ferromagnetism, such as the Hubbard 
model,— — there is no simple and comprehensive physi- 
cal picture for the mechanism that drives ferromagnetic 
order §r21 

The physical reason which hampers a straightforward 
understanding of itinerant ground-state ferromagnetism 
probably consists in the fact that the ordering and actu- 
ally the formation of local magnetic moments is a strong- 
coupling phenomenon and thus in general not capable by 
perturbative techniques. This is opposed to antiferro- 
magnetic order, for example>ii Slater or band antifer- 
romagnetism is accessible by weak-coupling approaches, 
Heisenberg or local-moment antiferromagnetism emerges 
in effective low-energy models. 

It is therefore important to recognize that the same 
problems already show up in the Hubbard model on 
infinite-dimensional lattices. This limit, however, is 
rigorously accessible by dynamical mean-field theory 
(DMFT)fi2r— and many insights concerning itinerant 
ferromagnetism could be obtained in this way.— First 
of all, quantum fluctuations are recognized as essen- 
tial. This means that a static mean-field approach, like 
Hartree-Fock theory, cannot grasp the main physics and 
largely overestimates the tendency to collective ordering. 
Further, ferromagnetism requires a strong local Coulomb 
interaction and must therefore be investigated by non- 
perturbative meansJ^"— In addition, however, subtle de- 
tails of the non-interacting electronic structure are like- 
wise important, e.g. a strong asymmetry of the local den- 
sity of statesi^"— Non-local parts of the Coulomb inter- 
action, a non-local ferromagnetic Heisenberg exchange 
coupling, for example, do affect the magnetic ground- 
state phase diagram but are found to be of lesser im- 



portance as compared to the Hubbard-C/ in general. Fi- 
nally, ferromagnetic order strongly competes with anti- 
ferromagnetism and is realized away from half-filling. 

While DMFT can be regarded as the optimal theoret- 
ical framework to deal with strong local quantum fluctu- 
ations and to understand their effect on itinerant ferro- 
magnctism, it is still a mean-field approach. This means 
that the feedback of non-local two-particle, e.g. magnetic, 
excitations on the one-particle spectrum and also on the 
thermodynamics is neglected. It is presently unclear, as 
how severe this approximation must be regarded when 
studying low-dimensional systems, for example. Spin- 
charge separation;^ to mention a prominent example of 
a non-local quantum effect in one-dimensional chains, 
cannot be described by DMFT. As concerns ferromag- 
netic order in one-dimensional itinerant systems, how- 
ever, there is reason to be more optimistic that DMFT 
may capture the essential physics: Namely, magnetic cor- 
relations and thus the feedback of magnetic correlations 
on the ground state can be expected to become less im- 
portant for fillings well below half-filling. At and around 
half-filling antiferromagnetic correlations dominate any- 
way. In fact, numerically exact studies by means of 
the density-matrix renormalization group (DMRG) for 
the ti-t2 one-dimensional Hubbard model^ir— yield a 
ground-state ferromagnetic phase diagram which shows 
striking similarities with the DMFT results and confirm 
the main qualitative results listed above. 

This situation has motivated the present study which 
employs the variational cluster approach (VCA)^i^ to 
investigate the ferromagnetic ground-state phase dia- 
gram of the infinite- and the one-dimensional Hubbard 
model. The VGA is a thermodynamically consistent^S 
cluster mean-field approach which determines the elec- 
tron self-energy by exploiting a general variational 
principle;^"— Different approximations can be con- 
structed by the choice of different reference systems that 
define the space of test self-energies for the variational 
principle. In this way, single-site mean-field approxima- 
tions, very close to the DMFT, as well as cluster approx- 
imations can be constructed which include the local but 
also non-local quantum fluctuations, respectively. 

By comparison with previous DMFT and DMRG re- 
sults it should be possible to answer the following interre- 
lated questions: How sensitive is a ferromagnetic state on 
a one-dimensional chain to local quantum fluctuations? 
What affects its stability more, local or short-range non- 
local fluctuations? Is a single-site mean-field approach 
sufficient to predict stable ferromagnetic phases? How 
much docs it improve compared to a purely static ap- 
proach? Docs an inclusion of short-range non-local fiuc- 
tuations improve the predictive power? Answers to these 
questions are particularly important for future studies of 
magnetic nanosystems in more complex geometries such 
as clusters, coupled chains, etc. and including more or- 
bitals per sites since those systems are in most cases 
not accessible to an exact numerical treatment via the 
DMRG. 



A second and likewise important goal of our study is 
to advance the variational cluster approach: Its evalua- 
tion requires the repeated calculation of the self-energy of 
the reference system for different one-particle parameters 
which serve as variational parameters. At zero temper- 
ature the numerical solution has to be performed using 
exact-diagonalization techniques. It is then clear that the 
quality of the approximation is limited by the exponen- 
tial growth of the reference system's Hilbert space, i.e. 
by the limited number of sites that can be taken into ac- 
count. However, an increasing number of sites in the ref- 
erence system at the same time means that the number of 
variational parameters increases. To study ferromagnetic 
phases, the number of parameters is doubled because of 
the additional spin-dependence of each of the parame- 
ters. The variationally determined ground-state energy 
becomes decreasingly sensitive with each additional vari- 
ational degree of freedom considered. This tightens the 
need for extremely accurate computations. Here our goal 
is to present and discuss different technical improvements 
of the VGA. 

The paper is organized as follows: In the following 
section |lT] we briefly review the theoretical concept and 
then address the different technical issues important for 
a reliable numerical evaluation in Sec. IIIII Results for 
the Hubbard model in inflnite dimensions and in one di- 
mension are presented and discussed in Sec. IIVI and a 
summary of the main conclusions is given in Sec. [V] 



II. VARIATIONAL CLUSTER TECHNIQUE 

We consider the single-band Hubbard modelS."— iu one 
dimension with nearest and next-nearest neighbor hop- 
ping ti and t2, respectively (except for Sec lIVBl) . Using 
standard notations, the Hamiltonian reads 



H = -t 
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where < • > and ^ • 3> restrict the independent sums 
over lattice sites i and j to nearest and next-nearest 
neighbors, respectively. <t —\, \. is the spin projection. 
The strength of the on-site Hubbard interaction is given 
by t/. In the following we consider finite Hubbard chains 
consisting of L sites and assume periodic boundary con- 
ditions. Unless stated differently, we set ii = 1 to fix the 
energy scale. 

Galculations are performed using the variational clus- 
ter approximation22i2^ (VGA) which is a quantum cluster 
mean-field approach based on the self-energy functional 
theory^-^ (SFT). Gentral to the SFT is the self-energy 
functional 



n[Yi\ = Trln(G, 



-1 





F\n 



(2) 



which provides an exact functional relation between the 
self-energy S (with elements Eijo-(a;ji)) and the grand 
potential Q, of Eq. ([T]) at temperature T and chemical 
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FIG. 1: Examples for reference systems with L^ correlated 
sites (blue circles) and one additional bath site (red circles) 
per correlated site generating single-site (mean-field) and clus- 
ter approximations. The (spin-dependent) on-site energies Ec 
and Eb as well as the (spin-dependent) hybridization V are 
optimized. Hopping parameters ii and t'2 are kept fixed at 
their physical values. 



potential ^. Furthermore, In denotes the main branch 
of the complex logarithm, and Tr = T^^ exp(ia;„0"'') tr 
with tr being the trace over the spatial degrees of freedom 
and 0^ a positive infinitesimal which ensures convergence 
of the sum over fermionic Matsubara frequencies a;„ . Go 
is the free single-particle Green's function of the model 
which can be assumed to be known. The functional 
F\Yi] is formally defined as the Legendre transform of the 
Luttingcr-Ward functional $[G] which in turn is defined 
diagrammaticall}*^ or through a functional integral.— 

The self-energy functional is constructed such that it 
becomes stationary at the exact (physical) self-energy of 
the model system Eq. ([1]): 



5VL[Ys] = , 



(3) 



Due to the fact that F[S] is not known explicitly, an ap- 
proximation must be employed to make use of this vari- 
ational principle. The idea of the SFT is to restrict the 
variation of the self-energy in Eq. ([3]) to a subspace of trial 
self-energies which is spanned by the exact self-energies 
of a certain reference system. On this subspace, the self- 
energy functional can be evaluated exactly, provided that 
the reference system has the same interaction part as the 
original model and provided that an exact (numerical) 
cornputation of the test self-energies is possible (see Refs. 
l30l - [33 for details). Gcnerically, a reference system is a fi- 
nite Hubbard cluster with the same Hubbard-C/ but with 
its one-particle parameters (hopping and on-site energies) 
serving to parametrize the test self-energies. Usually, one 
selects a limited number of one-particle parameters as 
variational parameters A — (Ai, . . . , A^). Then, the sta- 
tionarity condition Eq. ^ is approximated requiring the 
function 



n{\) = f][s(A)] 



(4) 



to be stationary, i.e. dQ,{\)/d\ — 0. 

For the present study of the one-dimensional Hubbard 
model, we consider chains with Lc < 5 correlated sites as 
reference systems. One additional uncorrelated ("bath") 
site is attached to each correlated site (i.e. Ug = 2 local 
degrees of freedom), see Fig. [TJ Physically, this means 



to take retardation effects into account, i.e. with this 
choice of the reference system, local (temporal) fluctu- 
ations are included to some degree. Local fluctuations 
are treated exactly in the limit of ris = 00 only, i.e. for a 
continuum of bath degrees of freedom (see Ref . l32i ) . This 
would correspond to the (cellular) dynamical mean-field 
(C-DMFT) approach^. It has been demonstrated in var- 
ious contexts)^"— however, that the main effect of local 
correlations is already accounted for with Ug ~ 2, i.e. 
the essential step is the one from a plain VGA (no bath 
sites) to an n^ = 2- VGA while more bath sites give sec- 
ondary corrections. This is important since only a finite 
(small) number of sites can be treated when using an 
exact-diagonalization technique at temperature T = 
to compute the chain self-energy. Physically, the main 
point is that the n^ = 2 reference systems already al- 
low for a local (Kondo-type) singlet formation to screen 
the local magnetic moments. A non-local singlet forma- 
tion is possible for reference systems with Lc > 2. This 
describes the feedback of non-local magnetic correlations 
on the single-particle excitation spectrum. The degree to 
which, quite generally, spatial correlations are accounted 
for by the VGA is controlled by the choice of Lc, ranging 
from a (dynamical) single-site mean-field approximation 
for Lc ~ I over cluster mean-field approaches to the exact 
solution that is (in principle) obtained with Lc = 00. 

For a given reference cluster (Fig.[T|), we treat the on- 
site energies of the correlated and of the bath sites, Sc and 
Sb, as variational parameters. This ensures thermody- 
namic consistencj^ with respect to the particle density 
n: Even within the approximation, exactly the same re- 
sult is obtained for n which either can be determined via 
a ^-derivative of the SFT grand potential at stationarity 
or via a frequency integral over the one-electron spectral 
density corresponding to G = {Gq^ — S)~^ with the op- 
timal self-energy. Furthermore, to control the temporal 
fluctuations, the hybridization V is considered as varia- 
tional parameter (see Fig. [T]). For the present study it 
is important to allow for a possible spin dependence of 
all variational parameters to have thermodynamical con- 
sistency with respect to the magnetization in addition. 
As a simplifying but excellent*^ assumption to limit the 
number of parameters. Sc, £b and V are taken to be site 
independent, and the hopping parameters within the ref- 
erence chain arc fixed at their original values, i.e. t'l = ti 
and t'2 = t2- 

Galculations are performed for the grand canonical en- 
semble keeping the chemical potential fj, fixed. Due to the 
discrete energy spectrum of the finite reference cluster, 
and due to the U(l) symmetry of the cluster Hamiltonian 
H', however, the cluster ground state reveals a fixed to- 
tal particle number TV' within finite ^ ranges. Therefore, 
the cluster electron density n' = N' /{2Lc) is a discontin- 
uous function of n which would give rise to discontinuous 
behavior of the self-energy and thus of all observables. 
This fact elucidates the second major motivation for in- 
troducing one bath site per correlated site: Bath sites 
serve as charge reservoirs. The entries i,j of the self- 



energy T,ija{io) are restricted to the correlated sites, and 
for a half-filled cluster, i.e. n' = 1 or A^' = 2Lc, the elec- 
tron density on the correlated sites can vary in the entire 
range from ?ij, = to n^ = 2 (see Ref. [3^ for a detailed 
discussion). This is important not only for studies of 
density dependencies but also for ferromagnetic phases. 
With the help of the bath sites, an arbitrary and continu- 
ous variation of the cluster magnetization m' 



''ct' 



'ci 



can be achieved in the same way. We will comment on 
this in the discussion of the results below. 



III. NUMERICAL EVALUATION 

The actual calculations are performed for finite Hub- 
bard chains (Hamiltonian H, Eq. ([T])) with typically 
L = 0(10'^) sites and assuming periodic boundary con- 
ditions. It is convenient to consider the reference system 
(Hamiltonian H') as being composed of Nk = L/Lc iden- 
tical and disconnected clusters consisting of Lc sites each, 
i.e. the sites of the reference system form a translation- 
ally invariant superlattice of Nk supersites. Via the self- 
energy, this superlattice structure is also imposed on the 
expectation values of observables of the original model. 

The grand potential at zero temperature can be calcu- 
lated asi^ 

f^(A) = o'-^c.;e(-c.:j+^c.„(fc)e(-c.„(fc)). (5) 



Here, il' is the grand potential of the reference system 
and uj'„ are the poles of the one-electron Green's func- 
tion of the reference system. These can be calculated 
exactly by means of the Lanczos approaches Note that 
the w^ do not depend on the wave "vectors" k of the first 
Brillouin zone of the superlattice as the clusters are dis- 
connected, and thus the fc-sum in the second term simply 
yields a factor Nk- w„(fc) are the poles of the (approxi- 
mate) Green's function of the model Eq. ([1]). Finally, Q 
denotes the Heaviside function, and the A dependence of 
fl' and of the poles w^ and uJnik) is implicit. 

There are several technical points which are essential 
for a reliable numerical evaluation of the VGA and there- 
with for the interpretation of the results. One of the ma- 
jor intentions of the present paper is to show how one can 
efficiently deal with the different finite-size effects in the 
evaluation of Eq. ([S]), in particular close to a second-order 
phase transition, and with the problem of finding station- 
ary points in a high-dimensional parameters space. 



A. Exact frequency summation 

The first problem consists in the infinite sums over 
Matsubara frequencies in Eq. ([5]). For not too large 
reference systems, this is performed conveniently and 
numerically exact by means of the so-called Q-matrix 
technique i^^i^^ Let A' be the diagonal matrix with the 



poles w^ of the cluster Green's function G' as diagonal 
elements. The Lehmann representation of G' can then 
be written as 



G'{uj) = Q 



1 



cj — A 



iQ' 



(6) 



with an appropriate weight matrix Q^ia),n- Note that 
QQ^ = 1 7^ Q^Q. Using this Lehmann representation in 
the definition Gk = (Gq ^ — S)~^ of the Green's function 
of the original model, it is easy to see that the poles uin{k) 
of Gk then can be obtained as eigenvalues of the matrix 



M(/fc)=A' + Qty(fc)g 



(7) 



Here V{k) = e{k) — t' is the difference between the one- 
particle parameters of the original and the reference sys- 
tem where the matrix e{k) is the Fourier transform of 
t with respect to the superlattice. In practice, the effi- 
ciency of the Q-matrix technique is set by the dimension 
of M{k), i.e. the number of poles of G', which, using 
Lanczos as a cluster solver, typically amounts to 0(100). 
For Lc < 8, the repeated diagonalization of M{k) for all 
k, and for Lc > 8 the Lanczos diagonalization of H' rep- 
resents the dominant contribution to the necessary total 
CPU time, respectively. 



B. Interpolative fc-summation 

The fc-summation is much more tedious. From Eq. ([5]) 
it can be read off that Q{X) is a non- analytic function for 
finite L. We first discuss the (approximate) one-electron 
excitation energies a;„(/c) of the original system. A sign 
change of one of the energies a;„(fc) as function of a vari- 
ational parameter A^ causes a kink in i7(A) due to the 
Q function. Such kinks have a negligible relative weight 
in the fc-suni and can be ignored in the thermodynamic 
limit i — >■ oo (if the interacting density of states at the 
Fermi edge stays finite). For finite L and in regions of 
the parameter space where n(A) is nearly fiat, however, 
the mentioned kinks may lead to artifacts or at least to 
severe convergence problems for numerical techniques to 
find stationary points, particularly if derivatives of i7(A) 
are required. 

In principle, this finite-size effect can be controlled by 
increasing the system size L and thereby the number of 
fc-points Nk- Although the computational effort is only 
linear in Nk, we found it to be much more effective to 
employ an interpolation algorithm which artificially in- 
creases the number of A:- vectors while keeping the system 
size fixed. For two adjacent fc-vectors fci and ^2, we inter- 
polate between the excitation energies a;„(fci) and a;„(fc2) 
instead of calculating a;„(fc) for intermediate k by diago- 
nalization of M{k) . Simple linear interpolation turns out 
to be sufficient. The effect is a smoothing of the function 
51(A) which considerably stabilizes the subsequent opti- 
mization procedure without a significant increase of the 
computational effort. 
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FIG. 2: Self-energy functional !^(S(F, £c, £&)) = ^{V, Ec, et), 
plotted versus the hybridization strength V while keeping the 
on-site energies fixed at their optimal values, Sc = —0.14 and 
£(, — —0.24. Calculation for ii = 1 (this sets the energy scale 
throughout the paper), i2 = 0, C/ = 4, p = —0.3 and using 
the Lc = 2 reference system. Solid lines represent calcula- 
tions without the interpolation method, dashed lines refer to 
calculations with a total (original plus interpolated) number 
of poles Nk,i„f 



Which pairs of poles at fci and ^2 correspond to 
each other, respectively, is actually unknown (as long as 
one does not analyze the corresponding eigenvectors of 
M{k)) but for practical purposes it is sufficient to sort 
the respective pole sets and assume the n-th pole at fci to 
correspond to the n-th pole at fc2, i.e. "level" crossing is 
excluded. If the number of fc-points is sufficiently large, 
the interpolation procedure affects contributions to the 
sum over k and n in the last term of Eq. ([S]) only in 
those cases where a;„(fc) crosses zero. Simple continuity 
arguments then show that a "level" crossing is unlikely in 
those cases, i.e. the possible error of disregarding "level" 
crossings is 0{1/Nk). This simple idea can also easily be 
generalized to higher dimensions. 

As an example to illustrate the interpolation scheme 
we show in Fig. [2] the self-energy functional as func- 
tion of the hybridization V for a reference system with 
Lc = 2. Convergence on the scale of the figure is ob- 
tained with Nk ~ 1000 (red solid line). As can be seen, a 
comparatively smooth curve can also be obtained with 
Nk = 100 original but additional 900 interpolated fc- 
points (dashed green line). The comparison shows that 
the trend oin{V) is essentially unaffected by the interpo- 
lation scheme. This merely produces a tiny overall shift 
of Q(V) which is irrelevant for the determination of min- 
ima and maxima. Using much less (yellow dashed) or 
no interpolated fc-points (blue solid line) introduces the 
above-mentioned kinks. For the example shown here, 
where the SFT grand potential is rather flat, these kinks 
would render a reliable determination of the optimal hy- 
bridization strength impossible. 



C. Level crossing in the reference cluster 

Let us now turn to the second possible source of a non- 
analytic behavior of r2(A), namely a sign change of one 
of the single-electron excitation energies uj'„ of the ref- 
erence system as a function of a variational parameter. 
Consider an electron-removal process, for example. Here 
w^ = Eo{N') - En{N' - 1) - M < where K(^') is 
the n-th excited eigenenergy in the invariant subspace 
of H' with (total) particle number N'. If, as a function 
of A, the excitation energy tu'^ — >■ 0, the ground state of 
the reference system becomes degenerate with an even- 
tually new ground state in the A'"' — 1 subspace. Hence, 
w^ = would indicate a level crossing and a discontinu- 
ous change of the ground state of the reference system. 
This in turn would induce a discontinuous change of the 
self-energy and thus a discontinuity of the SFT grand po- 
tential £7 (A) which was unphysical for obvious thermody- 
namical reasons. It is therefore of utmost importance to 
keep N' = const, i.e. to ensure that the stationary point 
of the SFT functional (and a finite environment in pa- 
rameter space) always corresponds to the same A^'. The 
same holds for z-componcnt of the total spin. 

In principle, one might try to ignore a sign change 
of w^j as a function of the optimal A, i.e. as a function 
of a model parameter, and formally calculate the self- 
energy from the ground state of a subspace with given 
N'. Besides the fundamental problem that this would 
actually correspond to a non-equilibrium situation, such 
a procedure also cannot work in practice: As the above 
discussion has shown, cj^^ = would then induce a strong 
kink in i7(A) which cannot be smoothed unless extremely 
large reference clusters with Lc — >■ oo are considered. 



D. Local optimization 

Tracing a stationary point of S7(A) as a function of a 
model parameter can be accomplished by a local tech- 
nique, i.e. assuming the stationary point Agt to be close 
to a starting point Aq. A naive application of New- 
ton's method to find a zero of Vr2(A) has turned out 
to be inefficient, however, since equipotential surfaces 
ri(A) = const, are usually highly anisotropic. Below 
we briefly describe our modified algorithm which uses 
an adaptive local coordinate frame with a directionally 
dependent calculation of partial derivatives. 

Let A„ denote the set of variational parameters at iter- 
ation step n. Assuming n{X) to be approximately given 
by a quadratic form, the next estimate is 



^n-l-l 



A„ - Jf-i [vn{x)]^__ 



where 



H„ 



d^n{\) 



dXidXi 



(8) 



(9) 



A=A„ 

is the Hessian at A„. A numerically stable evaluation of 
the Hessian (and the gradient) is crucial here. This can 



be achieved iteratively by principal axis transformation: 

(10) 



if„ = U^DJJ^ 



The diagonal matrix D„ contains the eigenvalues of if„. 
New coordinates A are defined via the orthogonal trans- 
formation 



A = C/iA 



(11) 



The Hessian for the next n + 1-st step is calculated in the 
new frame: 



Hn 



+ 1,13 
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dXidXj 



(12) 
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FIG. 3: Crossover from the Lc = 1 to the Lc = 4 reference 
system (left) and from the Lc = 2 reference system to the 
one with Lc = 4 (right). Dashed hnes represent intra-cluster 
hopping parameters scaled by a factor a. It is assumed that a 
stationary point of the SFT functional is known for the case 
where the intra-cluster hopping parameters are switched off 
(a = 0). The stationary point is traced locally while adiabat- 
ically switching on the hopping parameters, < a < 1. This 
yields a stationary point for the respective Lc > 1 reference 
system (a — 1). 



Here ri(A) = 0(A(A)). Inverse transformation 



Hn+l — UnH„ + lU 



(13) 



yields the new Hessian in the original frame which is re- 
quired for the next iteration step. The main point is that 
in principal coordinates we have 



n 



ix)-n{Xst) = ^J2 



d^n{\) 



dxl 



{K-\st.,f (14) 



A=A„ 



for A close to Ast- Hence, Hn+i becomes almost diago- 
nal and can be calculated as a difference quotient using 
discrete steps which depend on the principal direction: 



AA, 



An/H„ 



(15) 



Here AJ7 is a suitably chosen constant. This implies that 
partial derivatives along directions in parameter space 
where ri(A) is almost flat are computed with a large AAi, 
while a small AA^ is used along directions where 57(A) is 
strongly curved. This has turned out to be crucial for a 
numerically stable algorithm. 



E. Global optimization 

As a prerequisite for a local method to trace a station- 
ary point, a global method must be available which is ap- 
plicable even if a reasonable starting point is not known. 
Except for global minimization algorithms which can be 
applied to minimize |Vf2(A)p, for example, there is no 
general global technique to find stationary points in a 
multidimensional space unfortunately. 

In context of the SFT, however, there is an elegant 
solution to this problem since any local method can be 
converted into a global one with the help of a crossover 
procedure as has been pointed out in Ref. \^ The main 
idea is to modify the original system by switching off 
the intercluster hopping (in the same way as it is done 
in the reference system). For this truncated system, the 



VGA trivially yields the exact solution, and the station- 
ary point is trivially given by one-particle parameters of 
the reference system which are equal to those of the trun- 
cated one. One then adiabatically switches on again the 
inter-cluster hopping in the truncated system, i.e. one 
replaces iintcr -^ atintor and increases the parameter a 
from a = to a = 1. During this adiabatic process the 
stationary point of the reference system Ast(a) can be 
traced by means of a local optimization method. Finally, 
Ast = -^31(0^ = 1) is the stationary point of the original 
system. 

Here, we present a variant of this crossover trick which 
makes use of the fact that it has turned out to be easy to 
globally find a stationary solution for the Lc = 1 refer- 
ence system (i.e. for the dynamical impurity approxima- 
tion). An adiabatic crossover from the Lc = 1 reference 
system to an Lc > 1 reference system can be performed 
then by introducing a dimensionless parameter a to scale 
the nearest-neighbor and next-nearest-neighbor hopping 
in the reference system: t'l — >■ at'i and tj ^^ '^t2- For 
a = we recover the mean-field solution, for a = 1 
we have the VGA with Lc > 1. This procedure can be 
applied to cross over between two arbitrary reference sys- 
tems. Examples are given in Fig. 21 

For our calculations using the Lc = 4 reference sys- 
tem, we started from Lc = 1 and have changed a from 
to 1 in steps of 0.05. An example is shown in Fig. 
m The crossover procedure can be done along two dif- 
ferent routes, namely from Lc = 1 to Lc = 4 directly 
and, in two crossover steps, from Lc = 1 via Lc = 2 to 
Lc = 4. The resulting optimal variational parameters 
are the same for both routes as can be seen from the fig- 
ure. It is physically plausible that with increasing a and 
cluster size, the values of the optimal parameters tend to 
approach the "physical values" of the original system, i.e. 
V = and Sc = while for Lc = 1 stronger deviations 
are necessary to partially compensate for the effect of the 
truncated inter-cluster hopping. 
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F. Calculations at fixed density and magnetization 

Sclf-cncrgy-functional theory has originally been devel- 
oped using the grand-canonical ensemble. At zero tem- 
perature, starting from a grand-canonical Hamiltonian 
with chemical potential /i and magnetic field B, 



n = H 



A^ X^^^^t + "'^^ ^ ^ ^(''^■n - "-4) ' (16) 



the SFT grand potential H, = ri(A, /i,i?) is a function 
of the variational parameters A [see Eq. (j4])] and of ji 
and B (and other model parameters, such as U). The 
variational parameters are fixed by dVl{\,iJL,B)/d\ = 
while, at the respective stationary point, the derivatives 
with respect to ^ and B yield the expectation values of 
the total particle number and magnetic moment: 



{N) 



dn{X,n,B) 



m - - 



dn{X,n,B) 
dB 



(17) 



We also define the electron density n ~ N/L = 
^i„{niij) / L and the magnetization m = M/L = 
Y.ia Za{nia)/L with z-t-^j, = ±1. 

To construct phase diagrams it is much more conve- 
nient, however, to keep n instead of fi fixed, e.g. to study 
U dependencies at fixed density n. If there is no manifest 
particle- hole symmetry (off half- filling or for ^2 7^ 0), the 
corresponding chemical potential is not known a priori. 
Furthermore, it is highly desirable to perform calcula- 
tions at fixed m instead of B to search for a ferromag- 
netic phase: Starting from a paramagnetic solution with 
m = and adiabatically increasing to, one simply has to 
trace the solution and find the corresponding B = B{m). 
A spontaneous ferromagnetic solution is then indicated 
by a finite m with B{m) = 0. 



Consider the twofold Lcgendre transformation from 
the grand potential VI via the free energy F = D, + jiN to 
the Gibbs energy G = n + ^iN + BM: 



0(A, n, B) ^ F{X, N, B) ^ G(A, N, M) . 



(18) 



At given N and M, the SFT Gibbs energy G(A, A^, M) is 
obtained from G(A, iJ,,B,N,M) = ^(A, n, B)+^iN+BM 
via the original stationarity condition 



dG ^^ dn ,^ 

9A=0^9A=*^ 
and two additional conditions fixing fi and B 



dG 






(19) 

(20) 
(21) 



Hence, one simply has to consider /i and B in addition 
to A as variational parameters. 

To illustrate the method we have performed VGA cal- 
culations using the Lc = I reference system at fixed den- 
sity 7T, = 0.7 and different U close to a second-order phase 
transition, see Fig. [5] The magnetization is treated as 
a given quantity. The top panel shows the Gibbs free 
energy as function of to. For U = 19 this is a con- 
vex function as it is prescribed by thermodynamic sta- 
bility. Uc = 19.7 marks a critical point above which 
the Gibbs energy becomes thermodynamically unstable 
within a certain range of magnetizations. As can be 
seen in the top panel, the phase is locally unstable for 
—0.18 < m < 0.18 where the Gibbs energy is con- 
cave. A thermodynamically stable state is obtained via 
a Maxwell construction. This yields the solid line. Be- 
tween —0.32 < m < 0.32 the Gibbs energy is a con- 
stant which implies B = dG/dM = 0. Hence, an in- 
finitesimal field B = 0^ will produce a finite magne- 
tization m = 0.32. States with |to| < 0.32 can real- 
ized by macroscopy phase separation. We conclude that 
Uc = 19.7 marks a continuous transition from the param- 
agnetic state to a state with spontaneous ferromagnetic 
order. The function B{m) (lower panel in Fig. [5]) can 
be discussed analogously. Local thermodynamic insta- 
bility is indicated by a negative slope, and instability 
with respect to a Maxwell constructed state is indicated 
by the dashed line. Finally, the same physics can be seen 
by looking at the free energy (at T = equal to the 
ground-state energy E = (H) — BM) given as a func- 
tion of B (see inset for U ~ 21). Note that to can be 
computed as a derivative of G or via an integration of 
the spin-dependent local density of states which, due to 
the thermodynamical consistency of the SFT, yields the 
same result. 

It is interesting that the mean-field [Lc = 1) approach 
yields a critical interaction Uc = 19.7 which is rather 
close to the numerically exact result Uc — 18.5 obtained 
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FIG. 5: Gibbs free energy G (per site) and the external 
magnetic field B as functions of the magnetization va. Cal- 
culations have been performed at fixed density n = 0.7 for 
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Lc — 1 reference system {L — 4000). For [/ = 21 the dashed 
line shows the result of the actual calculation. This solution, 
however, is unstable for |77i| < 0.32, and a solution with lower 
G is obtained by a Maxwell construction (solid line). The 
inset shows the free energy F (per site) as function of B for 
U = 21. 



via density-matrix renornialization groupi^ Characteris- 
tic for a mean-field approach is the square-root behavior 
of the order parameter m{U) close to the critical point, 
m ex \/U — Uc- This can be seen in Fig. |6l where the 
magnetization is displayed as function of U for fixed den- 
sity n = 0.7. The inset shows a linear trend of m? close 

to Uc. 



G. Ferromagnetic susceptibility 

For the calculations of m{U), as displayed in Fig. [51 
the simultaneous optimization of 8 variational parame- 
ters is required, namely the spin-dependent on-site ener- 
gies Ec and Eb, the spin-dependent hybridization strength 
V plus fi and B. The number of parameters can be re- 
duced if one is interested in the phase boundaries only. 
Rather than tracing a spontaneously symmetry-broken 
solution, a second-order critical point can be found from 
the divergence of a suitably defined susceptibility. Here 
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FIG. 6: Magnetization m as a function of the interaction 
strength U. For U = 21 a, spontaneous magnetization mris = 
0.32 is obtained corresponding to the blue lines in Fig. [S] 
Inset: m^ as a function of U close to Uc- The red dashed 
line is a linear extrapolation to m = 0. Parameters of the 
calculation: see Fig. (5] 



we consider the homogeneous static magnetic suscepti- 
bility X = —(1/^) d'^F/dB'^. The most obvious way 
to calculate x is to apply a small external homogeneous 
magnetic field B and to look at the linear response m, 
i.e. X ~ limB-toTTi/B. For this case, it is convenient to 
consider B as fixed which implies that only 7 variational 
parameters (for Lc = 1 or Lc = 2) have to be taken into 
account. The result for the Lc = 1 reference system is 
shown in Fig. [7] (green lines). The divergence of x ^t 
Uc = 19.7 is consistent with the Uc extracted from the 
order parameter in Fig. [6] As it is typical for a mean- 
field approach x~^ (s^e inset) is a linear function of U 
close to Uc- The same holds for Lc = 2. Again this has 
to be expected as critical phenomena should not depend 
on the reference cluster size. 

For the calculation of the susceptibility, a variational 
optimization of spin-dependent variational parameters is 
actually not necessary as was recognized by Eder4^ This 
can be seen in the following way. Consider the free energy 
F = F{\, B) [Eq. p5)) ] where we suppress the TV depen- 
dence in the notation. Due to the stationarity conditions, 
dF{\, B)/d\ = 0, the optimal A can be considered as a 
function of B, i.e. A = A(i3). Therefore, 



d OF 



(A(B),B)=0 



dB d\ 
Carrying out the differentiation, we find 

d^F 



(22) 



^^^-(A(i?),5)-'^^(^) 



(A(B),B)=0. (23) 



d\d\^ "■ " ' dB dBdX' 
This is a linear set of equations which can be solved by 
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This calculation involves spin-dependent parameter opti- 
mization and to take the field i? as a variational param- 
eter in addition. 

Finite size effects play a crucial role for the calculation 
of the susceptibility in the critical regime. It turns out 
that the numerically most expensive calculation where 
the magnetization is kept fixed is most stable against 
finite size effects. On the scale of Fig. [7] converged results 
are obtained for a comparatively moderate system size of 
L = 4000 sites. 



IV. RESULTS AND DISCUSSION 



A. Ferromagnetism in the Hubbard model 



FIG. 7: Homogeneous static magnetic susceptibility x — 
dm/dB\B=o (per site) as a function of the interaction strength 
U using the Lc — 1 reference system (fi = 1, t2 = —0.2, 
n — 0.7). The inset shows x~^{U)- Results are obtained by 
three different techniques. Green lines: spin-dependent opti- 
mization of the variational parameters at fixed small field B. 
Blue lines: spin-independent parameter optimization accord- 
ing to Eq. (f26)l . Red lines: optimizing the Gibbs energy by 
varying the external field at fixed small magnetization. Con- 
vergence is obtained for different respective system sizes L as 
indicated. 



matrix inversion to get 

dXjB) _ 
dB 



d^F 



dXdX 



dBdX 



(24) 



Novir, the susceptibility 
-d^F{X{B),B)/dB^. Hence 



IS 



given by X = 



_d_ / dF{X{B),B)dX{B) 
dB V dX dB 



dF{X{B),B) 
dB 



(25) 

Using Eq. (|22p and the stationarity condition once more, 
we see that the first term does not contribute, and thus 



X 



d^F{X,B)dX{B) d^F{X,B) 



dXdB dB 



dB^ 



(26) 



where dX{B)/dB can be eliminated using Eq. (|24|) . Con- 
sequently, for the calculation of x it is sufficient to 
consider a paramagnetic state and to optimize spin- 
independent variational parameters only. This strongly 
reduces the computational effort. Once a paramagnetic 
stationary point is found, partial derivatives according 
to Eq. (|26|) and Eq. ((24)) have to calculated with spin- 
dependent parameters A in a final step. The resulting x 
as a function of U is shown in Fig. [7] as the blue lines. 

A third way to determine the susceptibility is to keep 
the magnetization m fixed at a small value and vary the 
field B to optimize the Gibbs energy G, see red line in Fig. 
[71 Here, a divergence of x is indicated by B{m 7^ 0) = 0. 



Applying the Hartree-Fock approximation to the 
single-band Hubbard model, one is lead to the Stoner 



criterionr 



UpoiO) > 1 , 



(27) 



for the existence of a ferromagnetic instability. There- 
with, the calculation of the free {U = 0) local density 
of states (DOS) po{uj) at the Fermi edge a; = can 
give first insights where ferromagnetism is likely to oc- 
cur. Conceptually, however, the Hartree-Fock approach 
is a static mean-field theory, and quantum ffuctuations 
are neglected altogether. If at all, reliable results can 
be derived for the extreme weak-coupling regime where 
ferromagnetism is unlikely to occur. 

Despite the simplicity of the Hubbard model, only a 
few rigorous results on ferromagnetism are available 4^"— 
The Mermin- Wagner theorem^i^ excludes spontaneous 
breaking of the SU(2) symmetry for finite tempera- 
tures and dimensions lower then three. For the one- 
dimensional case and nearest-neighbor hopping, Lieb and 
Mattis^ have shown that the ground state for any even 
number of electrons is always a non-magnetic singlet in- 
dependent of U. A ferromagnetic ground state is also 
excluded in the low-density limit n 1— >■ irrespective of 
U as has been argued by Kanamoriii His T-matrix ap- 
proach, however, must be based on the assumption that 
weak-coupling perturbation expansion converges. Lieb^ 
has shown that a ferromagnetic ground state is excluded 
for a bipartite lattice with nearest-neighbor hopping at 
half-filling n = 1 any U > independent of the dimen- 
sionality. As has been demonstrated by Nagaokai^ the 
ground state of the half-filled model with one hole added 
is fully polarized for C/ — 00 on bipartite lattices (for 
fee and hep lattices with negative hopping integrals) in 
three or higher dimensions. While criteria for the sta- 
bility of the fully polarized state for thermodynamically 
relevant dopings could not be obtainedr^"— it is possi- 
ble to reduce the parameter space left for a stable Na- 
gaoka state in the thermodynamic limit by different vari- 
ational approaches j^Si^ Mielke and Tasaki^"— proved 
the stability of ferromagnetism for special lattices, such 
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as the Kagomc lattice, for which there are dispersionless 
parts of the Bloch band ("flat-band ferromagnetism" ) . 
In these systems the Fermi sea is degenerate with ferro- 
magnetic states for [/ = 0, and ferromagnetism becomes 
stable for U > 0. Miiller-Hartmann^ has considered 
the one-dimensional model with a next-nearest-neighbor 
hopping such that the free band has two degenerate min- 
ima. In the low-density limit a metallic ferromagnetic 
ground state is obtained due to ferromagnetic exchange 
in a corresponding effective two-band model. Similarly, 
Tasaki^ constructed a one-dimensional Hubbard model 
with next-nearest neighbor hopping which has a ferro- 
magnetic (insulating) ground state at quarter filling and 
sufficiently strong U . 



B. Ferromagnetism in infinite dimensions 

A comprehensive but approximate approach to fer- 
romagnetic order is provided by dynamical mean-field 
theoryii2r2i In the past several studies have addressed 
the magnetic phase diagram of the Hubbard model on 
infinite-dimensional lattices where the DMFT becomes 
exact. At least two routes towards ferromagnetic or- 
der could be identified: (i) On a particle-hole symmet- 
ric hypercubic lattice ferromagnetism is realized for very 
strong Coulomb interaction U and fillings close to half- 
filling ii2r— This is reminiscent of the Nagaoka state j^ 
(ii) On the other hand, a moderate Hubbard-[/ is suffi- 
cient for lattices with a free DOS exhibiting a pronounced 



1.0 1 m — 1-^ 



asymmetry 



20-22 



Here a ferromagnetic ground state is ob- 



served in large regions of the U-n phase diagram. A 
simple mechanism for ferromagnetic order is not appar- 
ent although some understanding could be achieve d^' ^"i^^ 
by techniques and arguments related to the Hubbard- 
I approach.— The Stoner criterion turns out to be in- 
adequate. Furthermore, also a realization of flat-band 
fcrromagnctisnt^ can be found^ on a Bethe lattice with 
infinite coordination. 

We have performed calculations for different lattices, 
i.e. for different free DOS, respectively, using the self- 
energy-functional approach for a reference system with 
one correlated and one bath site only, i.e. Lc = 1, see 
Fig. [TJ This is referred to as the dynamical impurity ap- 
proximation (DIA) in the following. Let us recall that 
the DIA with an infinite number of bath sites would ex- 
actly correspond to DMFTi^ While local quantum fiuc- 
tuations are treated exactly within DMFT, the DIA is 
much simpler and includes some local fluctuations only. 
However, due to the presence of the bath site it allows 
for the formation of a local (Kondo-type) singlet. Our 
goal is here to test the DIA by comparing with available 
DMFT results for the ferromagnetic phase. This serves 
as a benchmark of the approximation. Furthermore, as 
a computationally cheap method, the DIA allows for a 
more comprehensive study of the phase diagram. 

We start with the flrst route (i) towards ferromag- 
netism and consider flUings close to half-fllling and very 
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FIG. 8: Polarization p = m/n (with m = n-|-— nj_ and n = n-f-f 
nj_) as a function of the doping 5 = 1 — n for the hypercubic 
lattice in infinite dimensions at (7 = 50. 5c is the critical 
doping. The energy scale is set by the variance a^ = 0.5 of the 
Gaussian free DOS. Results of the DIA (lines) are compared 
with data from full DMFT-NRG calculations (points) taken 
from Zitzler et al>i^ Dashed line: inverse homogeneous static 
susceptibility x~^ as obtained from the DIA. 



strong Coulomb interaction. Here we can compare with 
DMFT results by Zitzler et alJ^ which have been ob- 
tained using the numerical renormalization group (NRG) 
as an impurity solver. The calculations have been car- 
ried out for the hypercubic lattice in infinite dimensions. 
Using the same conventions as in Ref. [la, the free DOS 
is given by 






(28) 



The variance of the Gaussian DOS cr^ = 0.5 sets the 
energy scale. 

The results are displayed in Fig. [H For very strong U 
the DMRG-NRG data predict an almost fully polarized 
ferromagnetic state at low dopings 5 = 1 — n. Note that 
as a consequence of the tails of the free DOS, the ground 
state cannot be fully polarized in a strict sense.— How- 
ever, this exponentially small scale cannot be expected 
to be visible in the data. With increasing doping the 
system undergoes a continuous phase transition to the 
paramagnet at a critical doping 5c = 14.6%. The result 
of the DIA agrees well with the full DMFT and likewise 
predicts a continuous transition from a fully polarized 
state to the paramagnetic phase. On the rescaled plot in 
Fig. [HI the agreement is even quantitative. However, the 
critical doping for the phase transition {5c = 0.185) is 
significantly overestimated as compared to DMFT-NRG 
{5c = 0.146). As this means a stronger tendency towards 
ferromagnetism, one may conclude that the DIA under- 
estimates the effect of local quantum fluctuations. 

The DIA results are consistent in themselves: The 
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FIG. 9: Ground-state phase diagram U vs. filling n for the 
asymmetric free DOS given by Eq. (|3ip with half band width 
D = 2 and asymmetry parameter a = 0.98. The solid line 
refers to the DIA. DMFT-QMC results (points) are taken 
from Wahle et al.— 



magnetization m can be calculated via the spectral theo- 
rem from the spin-dependent one-electron spectral func- 
tion, or as the derivative of the optimal grand potential 
with respect to an external magnetic field. Both com- 
putations yield the same result as has been checked nu- 
merically and as is clear from the formalism)^ We also 
checked numerically that the Luttinger sum rule is ful- 
filled. Within the DIA this must be respected^ - in the 
paramagnetic but also in the ferromagnetic state. In the 
spin-polarized metallic phase there are two Fermi sur- 
faces with Fermi-surface volumes for a =f , | 



T^FS,. = ^e(M-e(k)-S,(0)) 



(29) 



where e(k) is the Bloch band dispersion and Scr(w) the 
fc-indepcndent self-energy. The Luttinger theorem then 
reads as 






duj pa-{uj) 



(30) 



For a local and real self-energy, the interacting local DOS 
can be written as Pcr{co) = pQ{iO + 11— Yig-{uj)), and 
the Luttinger sum rule reads p = poa- + ^a-iO). Here 
poa- is a (spin-dependent) chemical potential of the non- 
interacting system such that the spin-dependent particle 
numbers are the same as for the interacting system. 

Next we consider the second route (ii) towards fer- 
romagnetism and consider a moderate Hubbard-[/ but 
a free DOS with a pronounced asymmetry. Here we 
can compare with the results of Wahle et al.— who em- 
ployed the Hirsch-Fye quantum Monte-Carlo method as 
an impurity solver for DMFT. Calculations have been 
performed for finite but low temperatures and could be 
extrapolated to extract a ground-state U-n phase dia- 
gram which is shown in Fig. IH] (points). The Bloch band 
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dispersion e(fc) enters the DMFT (and also the DIA) 
via the free DOS only. Hence, instead of specifying the 
lattice structure and the hopping parameters, one can 
likewise start from a certain model free DOS as input for 
the DMFT calculation. This has the advantage that the 
effect of the asymmetr y o f the free DOS can be studied 
systematically. In Ref. [2l| the following model free DOS 
has been considered: 



po(w) = c 



/D2 



W^ 



D 



au! 



(31) 



Here, a is a parameter which controls the asymmetry 
while the variance stays constant. One can continuously 
tune the DOS from the symmetric case a = 0, corre- 
sponding to the semielliptic DOS of the Bethe lattice 
with infinite coordination, over an asymmetric DOS with 
more and more spectral weight peaked in the vicinity of 
the lower band edge, to a DOS with an inverse square- 
root divergence at the band edge for a = 1 eventually. 

(1 



Furthermore, in Eq. ^^, c = {1 + y/1 - a?)/{TTD) is 
a normalization constant, and D is the half band width 
which is set to D = 2 to fix the energy scale. The DMFT- 
QMC results in Fig. [9] correspond to a strongly asymmet- 
ric DOS characterized by a = 0.98. 

As is obvious from Fig.[SJ a ferromagnetic ground state 
is realized in large areas of the phase diagram. For low 
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fillings a moderate Hubbard-t/ is sufficient for ferromag- 
netism. With increasing n the phase boundary Uc(n) 
increases. Note that large U values cannot be accessed 
easily within the Hirsch-Fye QMC approach. It appears 
that this phase diagram is ruled by a mechanism that is 
completely different from the Nagaoka mechanism that 
has been suggested to rule the physics in case (i). Con- 
trary to the results shown in Fig. [51 ferromagnctism be- 
comes more likely with increasing doping 5 = 1 — n and 
persists down to very small fillings. 

Fig. ini also shows the result of our DIA calculation for 
a = 0.98 (solid line). Again, we find a convincing qualita- 
tive agreement with the full DMFT. The phase boundary 
Uc{n) shows the same trend but is systematically shifted 
towards higher interaction strengths. We attribute this 
difference partly to the very sensitive dependence of the 
results on the asymmetry parameter. This can be seen 
in Fig. [10] where the result for the ground-state phase 
diagram from DIA calculations for different asymmetry 
parameters a are given. It is obvious that for a close to 
unity a tiny change of a and thus of the free DOS results 
in a strong shift of the critical U. 

Using the DIA one can easily trace the evolution of the 
phase diagram as a function of the asymmetry parame- 
ter. As can be see from Fig. [TUl Uc{n) can be very small 
for a — > 1, i.e. for the case where the free DOS diverges 
at the lower band edge. For a < 1 the critical interaction 
becomes large and eventually Uc —>■ oo for n — )• 0. With 
increasing n, however, the phase boundary soon devel- 
ops a minimum at n„iin and then becomes an increasing 
function of n. This minimum is located at low fillings 
for asymmetry parameters close to unity but then shifts 
to higher fillings for a less asymmetric free DOS. At the 
same time, C/c('^min) increases strongly. For a = 0.5 we 
find ?7inin ~ 0.9, and the ferromagnetic phase is confined 
to a small filling range close to half-filling and very strong 
Coulomb interaction. 

It appears that the two routes towards ferromagnetism 
(Nagaoka vs. asymmetry of the free DOS) are linked 
continuously. For even smaller asymmetry parameters 
a < 0.5 ferromagnetism disappears completely. The sym- 
metric case a = corresponds to a Bethe lattice with 
infinite coordination with a symmetric free DOS. Here a 
ferromagnetic state cannot be stabilized. This is again 
consistent with full DMFT (NRG) calculations,^ Obvi- 
ously, the stability of the ferromagnetic state not only 
depends on the asymmetry a but is also strongly affected 
by the detailed form of the symmetric free DOS since, 
as has been discussed above, for the symmetric Gaus- 
sian free DOS corresponding to the hypercubic lattice, 
there is again a ferromagnetic phase close to half-filling. 
It is an open question whether the latter can really be 
attributed to the Nagaoka mechanism. On the one hand, 
the Nagaoka mechanism needs closed loops on the lat- 
tice which are present for the hypercubic one but absent 
for the Bethe lattice. On the other hand, the DMFT is 
sensitive to the lattice structure via the free DOS only. 

Concluding, we find that the DIA gives qualitatively 



reliable results for the ferromagnetic ground-state phase 
diagram in all cases studied, as has been corroborated by 
the comparison with different full DMFT calculations. A 
reference system with a single bath site appears to be 
sufficient to capture the main physics although quanti- 
tatively there is a tendency to overestimate the range 
where ferromagnetism is possible. While local quantum 
fluctuations are included in the DIA in a very simple 
way only, the approximation allows for the formation of 
a local (Kondo-type) singlet. Together with the internal 
thermodynamical consistency of the approach and with 
the fact that Luttinger's sum rule is respected, this en- 
sures a reliable mean-field description of ferromagnetism. 



C. Ferromagnetism in one-dimensional chains 

Ground-state ferromagnetism in the one-dimensional 
Hubbard model is restricted by the Lieb-Mattis 
theorem^ which excludes a finite order parameter in case 
of nearest-neighbor hopping only, irrespective of the in- 
teraction U. Including a ncxt-nearest-neighbor hopping 
t2, however, ferromagnetism is proven to exist for U — oo 
in the limit t2 ^ (t2 < 0) for all densities-^^i^ This 
limit has to be contrasted to the limit ti = 0, but finite 
t2 (two-chain model) where the Lieb-Mattis theorem ap- 
plies again. In the low density limit, the ground state 
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FIG. 11; Magnetic ground-state phase diagram U vs. t2 of the 
one-dimensional Hubbard model for quarter fiUing (n = 0.5) 
as obtained by the dynamical impurity approximation (DIA, 
blue). The shaded area represents the parameter regime for 
which a continuous transition is found, fi = 1 fixes the en- 
ergy unit. Only the range t2 < is considered. Note the 
non-linea r scales w ith — oo < i2 < and < U < oo. 
Wcs = \/2tj + 2tl is the effective band width of the free DOS. 
The orange line ("Nagaoka") marks the critical interaction 
strength below which the fully polarized state becomes insta- 
ble against the paramagnet. The green line ( "HF/Stoner" ) is 
the phase boundary according to the Stoner criterion. 
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is ferromagnetic for t2 < —1/4 at U — oo^ With a fi- 
nite next-nearest-neighbor hopping, ferromagnetism oc- 
curs in a rather large part of the U-n phase diagram as 
has been demonstrated by DMRG calculations of Daul 
and Noack,24i25 

To study the effect of local and of short-range non-local 
quantum fluctuations on the stability of the ferromag- 
netic ground state and to test the predictive power of 
mean-field and cluster mean-field approaches, we have 
applied the dynamical impurity approximation (DIA) 
and the variational cluster approximation (VGA) (see 
Fig. [T]) to the model with t2 ^ 0. Note that a finite i2 
translates into an asymmetric free DOS and that ^2 7^ 
implies magnetic frustration with respect to antiferro- 
magnetic order. We exclusively consider ^2 < (ti = 1 
sets the energy scale) which implies that ferromagnetic 
order is expected to show up for fillings below half-filling. 

We first check whether or not the Lieb-Mattis theorem 
is respected by the most simple DIA. To this end, DIA 
calculations have been performed to map out the U-t2 
phase diagram at a fixed filling n = 0.5 (quarter filling). 
Calculations are done using chains with up to 8000 sites. 
The critical interaction for ferromagnetic order Uc is de- 
termined by the divergence of the paramagnetic suscepti- 
bility. X is calculated by using a finite but small external 
field {B = 0.01). The values for Uc obtained in this way 
are checked for selected <2 by calculating the external 
field for a given small magnetization {m = 0.01). Devi- 
ations are small, i.e. invisible on the scale of the figures 
discussed below, and can be neglected. 

The resulting phase diagram is shown in Fig. [TTJ The 
DIA predicts a Uc which varies strongly with \t2\- To be 
able to display the results for < [/ < 00 and — 00 < 
^2 < in a single picture, non-linear scales for ^2 and 
U have be en used. The effective bandwidth defined as 
Weff = ■\/2t^ + 2t|, i.e. the standard deviation of the 
free DOS, and ti are chosen as the relevant scales for U 
and ^2- 

In both cases where the Lieb-Mattis theorem holds, 
for t2 = and ^2 — > 00 the DIA predicts Uc —?' 00, i.e. 
the absence of ferromagnetic order. On the other hand, 
the static mean-field theory is clearly at variance with 
the exact theorem as can be seen from Fig. [11] where 
the green line ( "HF/Stoner" ) displays the divergence of 
the Hartree-Fock susceptibility as determined from the 
Stoner criterion. Note that the discrepancy between the 
DIA and the HF results is drastic except in the vicinity 
oft2 ~ -0.3 (|i2|/(|ii| + |i2|) «0.23) where the chemical 
potential of the non-interacting system coincides with a 
van Hove singularity in the free DOS. Here the Stoner 
criterion correctly predicts the ferromagnetic instability 
of the ground state. 

As one cannot expect that the Lieb-Mattis theorem 
is respected rigorously within a mean-field approach, we 
also performed calculations for different fillings. In fact, 
for n = 0.4 a divergence of the susceptibility is found for 
^2 = 0. However, the large critical value for the interac- 
tion, Uc ~ 30, indicates that the violation of the Lieb- 
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FIG. 12: Maxwell construction for the discontinuous phase 
transition for n = 0.5 and i2 = -0.34 (|i2|/(|fi| + |f2|) ~ 0.25) 
at the critical interaction Uc = 0.43. The shaded areas have 
the same size, i.e. AG = (see Eq. (|32|) '). There is a solu- 
tion with spontaneous ferromagnetic order at m, = 0.20. The 
second ferromagnetic solution at m f» 0.14 shows a negative 
susceptibility x = dm/dB < and is thus locally (and glob- 
ally) unstable. 
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FIG. 13: Gibbs free energies G of the paramagnetic and of 
the ferromagnetic solutions (top panel) and the correspond- 
ing magnetizations m (bottom) as functions of the interaction 
strength U at t2 — —0.34 and quarter filling. The red line 
shows the actual course of the stable solution. The magne- 
tization vanishes continuously at U — 0.49, and at the same 
interaction strength the susceptibility x diverges. The true 
phase transition is discontinuous and takes place at Uc = 0.43. 
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FIG. 15: The same as Fig. ll4l but using the variational cluster 
approach with reference systems shown in Fig. [T] with 1 < 
I/c < 4 (and the same number of bath sites, i.e. ris = 2) in 
the range 0.7 < n < 0.8. 



Mattis theorem is "weak" in the sense that it occurs for 
extremely strong interactions only. 

The comparison with Hartree-Fock theory tells us that 
local quantum fluctuations are very important. On the 
other hand, the comparison with VGA results shows that 
non-local fluctuations are not important in first place. 
For example, for ^2 = and for a very strong interaction, 
U = 10*, we find a divergence of the susceptibility at a 
critical filling n = 0.49 within the DIA while n = 0.42 
within the VGA for Lc ~ 2. This is the correct trend 
as the critical filling must vanish for Lc ^> oo due to 
the Lieb-Mattis theorem. As compared to the improve- 
ment of the DIA with respect to static mean-field theory, 
however, this appears as marginal. 

In most cases we find the phase transition to be discon- 
tinuous. Fig. ll2| gives an example. Here the homogeneous 
magnetic field B is shown as a function of the magnetiza- 
tion m for t-z = —0.34 and U = 0.43. This corresponds to 
\t2\/{\ti\ + \t2\) ~ 0.25 and U/{U + Wee) « 0.22. Sponta- 
neous ferromagnetism requires a finite order parameter 
TO at _B = 0. There are three solutions: (i) the param- 
agnetic state at to = which shows a positive suscep- 
tibility X = dm/dB, (ii) a thermodynamically unstable 
ferromagnetic solution with negative x, and (iii) a stable 
ferromagnetic solution with x > and m = 0.20. The 
point |t2|/(|ti| + |<2|) « 0.25 and U/{U + Wcs) ~ 0.22 
is just the transition point as can be seen from the area 
under the B(m) curve, i.e. from the Maxwell construc- 
tion, but is somewhat below the point which is plotted 
in Fig. [TT] and at which the susceptibility x diverges: 
Uc/{Uc + W,«) w 0.24. 

This explains itself in Fig. [T3l where the Gibbs free en- 
ergies and the magnetizations of the three solutions are 
shown for fixed ^2 as a function of U. Note that for T = 
and B = the Gibbs free energy is the ground-state en- 



ergy: G ^ E. It can be seen that the thermodynamically 
unstable solution always has the highest Gibbs free en- 
ergy. The actual phase transition therefore takes place at 
the interaction strength U = 0.43 where the Gibbs free 
energies of the (stable) ferromagnetic and of the para- 
magnetic solutions are crossing. This is consistent with 
the Maxwell construction in Fig. [12] since the difference 
in the Gibbs free energies of the paramagnet and the fer- 
romagnet is given by 
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On the phase boundary AG — 0. The divergence of 
X, however, is related to the continuous vanishing of the 
order parameter of the unstable solution at Uc = 0.49 (see 
Fig. [131 bottom). Note, that on the scale used in Fig. 
[m the difference between Uc and the true (first-order) 
transition points is almost invisible. 

The t2-range in which the transition is continuous is 
marked as the shaded area in Fig. [TTJ In addition, the 
figure shows the ^2 dependence of the interaction strength 
at which the fully polarized ("Nagaoka") state becomes 
unstable as compared to the paramagnetic state. This 
line crosses the phase transition line (diverging x) at 
|i2| ~ 0.10 (corresponding to |t2|/(|ii| + 1*21) « 0.09) and 
|i2| ~ 0.39 (corresponding to |i2|/(|ii| + 1*21) ~ 0.28). 
This implies that for jtal > 0.39 and for 1^2 1 < 0.09 the 
first-order transition is a transition from the paramag- 
netic to the fully polarized ferromagnetic state while in all 
other cases the magnetization jumps to a non-saturated 
value at the respective transition point. Our DIA results 
are consistent with the DMRG calculations of Daul^^ 
which yield a second-order transition at ^2 = —0.2 and a 
first order transition at t2 = —0.8 for quarter filling. 

For a systematic comparison of the results of the dy- 
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namical impurity approach with DMRG data?^ we fix 
the next-nearest-neighbor hopping to ^2 = —0.2 and map 
out the phase diagram U vs. fiUing n. The result is shown 
in Fig. [TJ] in comparison with static mean-field theory. 
The phase diagram turns out to be qualitatively similar 
to the result for infinite dimensions (see Fig.|9]). The crit- 
ical interaction Uc strongly varies with n and becomes ex- 
tremely large for fillings close to half-filling. Despite the 
simplicity of the reference system, the agreement with 
the DMRG date is reasonable. 

Within static mean-field theory, local magnetic mo- 
ments are formed in the ferromagnetic state only. This 
might be the right picture for the low-density regime. For 
fillings n > 0.5, however, static mean-field theory fails to 
reproduce the phase diagram as the tendency towards 
ferromagnetic order is overestimated drastically. 

Opposed to static mean-field theory, the DIA allows 
for local-moment formation already in the paramagnetic 
state and captures the correlated mean-field physics of 
the paramagnetic Mott transition at half-fillingi^i^li^ 
The high values for the Hubbard interaction necessary 
to produce ferromagnetic order can be understood in 
a picture where the ferromagnetic state evolves from a 
highly correlated paramagnet with preformed but disor- 
dered local magnetic moments. With increasing U and 
with increasing fillings n the local magnetic moments are 
more and more efficiently screened by a collective (single- 
band) Kondo effect. The tendency to screen the local mo- 
ments counteracts the formation of magnetic order and 
thus results in very strong critical interactions. Such a 
mechanism is already included in the DIA. From the rea- 
sonable agreement with the DMRG data and the strong 
improvement with respect to static mean-field theory, we 
therefore infer that this mechanism is essential. While 
the correct Kondo scale cannot be captured with a sin- 
gle bath site (ug = 2), the possibility to form a local 
singlet within a thermodynamically consistent approxi- 
mation appears to be a key ingredient to understand the 
phase diagram. 

Besides a screening of the local moment by local fluc- 
tuations, a screening by non-local fluctuations is conceiv- 
able. This would lead to non-local singlets - or even to 
long-range anii-ferromagnetic order. Consequently, such 
a mechanism is expected to be effective for fillings close 
to half-filling where non-local antiferromagnetic correla- 
tions are important. Note, however, that due to the Lieb- 
Mattis theorem the necessity to include a finite ^2 already 
suppresses antiferromagnetic order by magnetic frustra- 
tion to some degree. This might explain that the effect of 
non-local fiuctuations appears to be comparatively weak 
for intermediate fillings and significant for fillings close 
to half-filling only. 

This can be seen in Fig. [15] where we compare the 
DIA phase diagram with the results obtained from VGA 
calculations with finite clusters as reference systems: 
Lc ~ 2 — A while the description of the local degrees 
is unchanged {us = 2, one bath site per correlated site). 
For n < 0.75 the critical interaction does not change 



much while for n = 0.8, the critical U is strongly reduced 
in the cluster approach. The VGA thereby improves the 
agreement with the DMRG. 

Although the step from Lc = I to Lc > 1 appears to 
be essential close to half-filling, the results of the clus- 
ter approach have to be interpreted with some care since 
the expected convergence with increasing cluster size can 
hardly be seen for Lc < 4. This reflects finite-size er- 
rors the size of which can be estimated by comparing the 
results for different Lc among each other. Within this 
(considerable) error there is agreement with the DMRG 
results. On the other hand the VGA, and more impor- 
tant even the DIA, is able to predict the qualitatively 
correct trend for the phase diagram. It is also important 
to note within the DIA it is much easier to find, stabi- 
lize and trace magnetic solutions. For the clusters with 
ic > 1 we have not been able to find solutions in the 
entire filling range for reasons discussed in Sec. IIII Gl 



V. CONCLUSIONS 

The self-energy-functional theory has been applied to 
the Hubbard model in infinite and in one dimension to 
investigate spontaneous ferromagnetic order. Using dif- 
ferent reference systems generating different single-site 
and cluster mean-field approximations, i.e. dynamical im- 
purity and the variational cluster approximations, it is 
possible to study the effects of local and of short-range 
non-local quantum fluctuations on the stability of the 
ferromagnetic ground state. 

We find that local fluctuations are of crucial impor- 
tance to get a qualitatively correct phase diagram and 
to respect the Lieb-Mattis theorem. Opposed to static 
mean-field theory, ferromagnetic order quite generally 
requires substantially higher interaction strengths and 
can be understood as evolving from preformed but disor- 
dered local magnetic moments. The extremely large crit- 
ical interactions found in one-dimensional chains could 
then be attributed to the screening of the local moments 
by local (Kondo-type) correlations which becomes more 
and more effective for increasing filling or interaction 
strength. This local singlet formation, on a qualitative 
level, is already included in the DIA. Singlet formation 
due to non-local correlations, included in the VGA, ap- 
pears to be relevant for fillings close to half-filling only, 
while a ferromagnetic ground state can be obtained in 
large areas of the parameter space and down to the low- 
density limit in particular. The limited importance of 
(antiferromagnetic) non-local correlations is of course in- 
terrelated with the frustration of antiferromagnetic or- 
der due to a next-nearest-ncighbor hopping i2, with the 
Lieb-Mattis theorem, and, in the case of infinite dimen- 
sions, with an asymmetric free DOS. There is an ob- 
vious similarity of the magnetic phase diagram of the 
one-dimensional model with the phase diagram in infinite 
dimensions which again suggests that local correlations 
play the predominant role for ferromagnetic order. 
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This is of some importance for future investigations 
of ferromagnetism in nano-sized objects, e.g. chains or 
clusters on non-magnetic substrates, as it opens a route 
to study those systems by (dynamical) mean-field meth- 
ods which, as concerns the system geometry, are more 
flexible than the density-matrix renormalization group, 
for example. The situation may be contrasted, e.g., with 
the absence of long-range antiferromagnctic order in one- 
dimensional systems which is caused by non-local quan- 
tum fluctuations. In the latter case any mean-field ap- 
proach would be questionable a priori. 

The study of infinite-dimensional lattices using the 
DIA has shown that the previously known parameter 
ranges that are favorable for ferromagnetic order, namely 
low to intermediate fillings and moderate U in case of a 
strongly asymmetric free DOS and fillings close to half- 
filling and extremely strong U in case of a symmetric free 
DOS, are linked continuously. This demonstrates the dif- 
ficulty to find simple "mechanisms" for ferromagnetic or- 
der in the Hubbard model. Ferromagnetism should there- 
fore be seen as a complex phenomenon the description of 
which necessarily requires non-pcrturbative and thcrmo- 
dynamically consistent many-body techniques. 

To study ferromagnetism within self-energy-functional 
theory, at least six variational parameters have to be op- 
timized simultaneously, namely the spin-dependent one- 
particle energies of the correlated and of the uncorre- 
lated sites and the spin-dependent hybridization strength 
in addition. Two more variational parameters must be 
considered for calculations at fixed filling and magneti- 
zation which is convenient for the construction of phase 
diagrams. This can be accomplished by a number of 
technical improvements concerning (i) an accurate treat- 



ment of k- and frequency summations, (ii) optimization 
algorithms which adapt to the local structure of the func- 
tional, (iii) global optimization algorithms to find a sta- 
tionary point of the functional. For the calculation of the 
static and homogeneous paramagnetic susceptibility, an 
optimization of spin- independent parameters is sufficient. 

The comparison with dynamical mean-field theory and 
with density-matrix renormalization-group calculations 
for the infinite-dimensional and for the one-dimensional 
model, respectively, has been essential to rate the ap- 
proximations. For both, infinite dimensions and one di- 
mension, a simple DIA turns out to be sufficient for a 
qualitative and rough scan of the phase diagram. This 
might be sufficient in view of the fact that the Hubbard 
and similar models themselves represent strong simplifi- 
cations as compared to a real material. In one dimension, 
a cluster approach including short-range correlations, i.e. 
the VGA, appears to be necessary for fillings close to half- 
filling. A satisfactory convergence with increasing cluster 
size, however, could not be obtained. For future stud- 
ies of more complicated low-dimensional geometries, we 
therefore suggest to use the DIA in those ranges of the 
parameter space where there are no significant deviations 
from results obtained by the cluster approach. 
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